program main
  
  implicit none
  integer :: i,j
  integer, parameter  :: n = 6
  integer, parameter  :: mode_1 = 1
  integer, parameter  :: mode_2 = 1
  double precision :: omega_arr(4,5),B_imag_arr(4,5)
  double precision :: x,y,t
  double precision :: du,dv,dut,dvt
  double precision :: mu,c,bi
  
  omega_arr(1,1:5) = (/6.712844035888430d0,10.598136345660269d0,13.113329869474558d0,16.021054935234528d0,17.730944523371178d0/)
  omega_arr(2,1:5) = (/2.174059402269135d0,3.475669329435217d0,4.373718805991318d0,5.375869130532646d0,6.404312117277238d0/)
  omega_arr(3,1:5) = (/0.690252461874728d0,1.104283873555503d0,1.392820917483943d0,1.707902134335568d0,2.031628962752277d0/)
  omega_arr(4,1:5) = (/0.218370110186628d0,0.349373578609341d0,0.440766027573439d0,0.540340250707020d0,0.642652150833394d0/)
  
  B_imag_arr(1,1:5) = (/0.259859421131935d0,0.551757806713995d0,-1.111726166772277d0,0.958432243767886d0,-0.732513991983086d0/)
  B_imag_arr(2,1:5) = (/0.008979356435568d0,0.061797475082992d0,-0.195334916685379d0,0.716217080313942d0,-1.855812650827352d0/)
  B_imag_arr(3,1:5) = (/1.800969171774775d-5,1.608778558352172d-4,-5.996559491070863d-4,0.002791241585515d0,-0.010369491587158d0/)
  B_imag_arr(4,1:5) = (/1.946610327968597d-8,1.791032555369920d-7,-6.806632052120861d-7,3.250907109941674d-6,-1.252501580039413d-5/)
  
  bi = B_imag_arr(mode_1,mode_2)
  c =  omega_arr(mode_1,mode_2)/dble(n)

  t = 1.0d0
  x = 0.1d0
  y = 0.1d0
  mu = 1.d0

  call surface_wave_convex(du,dv,dut,dvt,x,y,t,mu,c,bi)
  write(*,*) 'If correct, this should return: 4.9261646866547626E-005  -1.3466800429152382E-004'
  write(*,*) x,y,du,dv
  
end program main
